﻿using LxBasic;

namespace LxCore
{
    public class TurbulencekOmega : ZoneHead
    {
        //常数来自wilcox论文
        public static readonly double Clim = 7 / 8.0;
        public static readonly double alpha = 13 / 25.0;
        public static readonly double betastar = 9 / 100.0;
        public static readonly double two_div_betastar = 2 / betastar;
        public static readonly double sigma = 1 / 2.0;
        public static readonly double sigmastar = 3 / 5.0;
        public static readonly double sigmado = 1 / 8.0;
        public static readonly double betao = 0.0708;
        public static readonly double PrT = 8 / 9.0;
        public static readonly double one_div_PrT;
        //红皮书
        public static readonly double Cmu14 = Math.Sqrt(Math.Sqrt(0.09));
        public static readonly double Karman = 0.41;
        public static readonly double Karman_Reciprocal = 1.0 / Karman;
        public static readonly double B = 5.25;
        public static readonly double Cbeta1 = 0.075;


        public NSEquation k;
        protected double[] k_x;
        public NSEquation w;
        protected double[] w_x;
        protected double[] miu_T;
        protected TurbBoundTypes TurbBoundType;
        protected CalBoundBase BoundValueCal;
        //share
        public ArrayScalar.BaseFace TotalViscosity;
        protected SIMPLEVelocity VelVector;
        protected Vector3d[] U_grad;
        protected Vector3d[] V_grad;
        protected Vector3d[] W_grad;
        protected IField<double> Density;
        protected IField<double> DynamicViscosity;
        protected WallDistance WallDistance;
        public override void Reuse(ZoneHead oldOne)
        {
            base.Reuse(oldOne);
            if (oldOne is TurbulencekOmega kw)
            {
                k = kw.k;
                w = kw.w;
            }
        }

        public override void CheckCreate(FvZone zone)
        {
            throw new NotImplementedException();
            //base.CheckCreate(zone);
            //LxBasic.CheckCreator.Instance(ref k, out NSEquation kns);
            //k.CheckCreate(zone);
            //this.k_x = k.XInn;
            //LxBasic.CheckCreator.Instance(ref w, out NSEquation wns);
            //w.CheckCreate(zone);
            //this.w_x = w.XInn;
            //LxBasic.CheckCreator.ArrayStrict(ref miu_T, zone.Cells.Length);
        }
        public override bool IsInitialized => k.IsInitialized && w.IsInitialized;
        public override void Initialize()
        {
            if(IsInitialized) return;
            if (!k.IsInitialized)
            {
                k.Initialize();
                k.LimitScalarValue();
            }
            if (!w.IsInitialized)
            {
                w.Initialize();
                w.LimitScalarValue();
            }
        }
        public void SetZone_Main(SIMPLEVelocity velocityVector, ArrayScalar.BaseFace gamma, IField<double> density, WallDistance WallDis, TurbBoundTypes turbBoundType)
        {
            this.TotalViscosity = gamma;

            this.VelVector = velocityVector;
            this.U_grad = velocityVector.U.Grad;
            this.V_grad = velocityVector.V.Grad;
            this.W_grad = velocityVector.W.Grad;

            this.Density = density;
            this.DynamicViscosity = velocityVector.DynamicViscosity;

            this.WallDistance = WallDis;
            this.TurbBoundType = turbBoundType;
            switch (turbBoundType)
            {
                case TurbBoundTypes.LowYPlus: BoundValueCal = new LowYPlusBound(); break;
                case TurbBoundTypes.HighYPlus: BoundValueCal = new HighYPlusBound(); break;
                case TurbBoundTypes.AllYPlus: BoundValueCal = new EnhanceYPlusBound(); break;
            }
            BoundValueCal.SetTurb(this);
        }
        //public void SetZone_part2(ValueBase<double> k_S, IFieldHost<double> k_iniValue, IFieldHost<double> k_relax, ConvectionSchemes k_convect, double k_minLim, ValueBase<double> w_S, IFieldHost<double> w_iniValue, IFieldHost<double> w_relax, ConvectionSchemes w_convect, double w_minLim)
        //{
        //    k.SetZone_CommonTerm(k_S, k_iniValue, k_relax);
        //    k.SetZone_DiffTerm(TotalViscosity);
        //    k.SetZone_ConvTerm(ConstField.Double1, k_convect, VelVector);
        //    k.SetLimitRule(true, true, k_minLim, false, double.MaxValue);
        //    w.SetZone_CommonTerm(w_S, w_iniValue, w_relax);
        //    w.SetZone_DiffTerm(TotalViscosity);
        //    w.SetZone_ConvTerm(ConstField.Double1, w_convect, VelVector);
        //    w.SetLimitRule(true, true, w_minLim, false, double.MaxValue);
        //}
        public void SetZone_kwEquation_ByTurbIntenAndVisRate(
            IField<double> k_S_b, IField<double> k_S_a, IField<double> iniTurbInten, IField<double> iniVelRef, IField<double> k_relax, ConvectionSchemes k_convect, double k_minLim,
            IField<double> w_S_b, IField<double> w_S_a, IField<double> iniVisRate, IField<double> w_relax, ConvectionSchemes w_convect, double w_minLim)
        {
            throw new NotImplementedException();
            //BoundkValue boundk = new BoundkValue();
            //boundk.SetParameter(iniTurbInten, iniVelRef);
            //k.SetZone_CommonTerm(k_S_b, k_S_a, boundk, k_relax);
            //k.SetZone_DiffTerm(TotalViscosity);
            //k.SetZone_ConvTerm(ConstField.Double1, k_convect, VelVector);
            //k.SetLimitRule(true, true, k_minLim, false, double.MaxValue);
            //BoundwValue boundw = new BoundwValue();
            //boundw.SetParameter(iniVisRate, k_x, Density, DynamicViscosity);
            //w.SetZone_CommonTerm(w_S_b, w_S_a, boundw, w_relax);
            //w.SetZone_DiffTerm(TotalViscosity);
            //k.SetZone_ConvTerm(ConstField.Double1, w_convect, VelVector);
            //w.SetLimitRule(true, true, w_minLim, false, double.MaxValue);
        }
        public void SetBound_WallOrSymmetry(FvBound thread)
        {
            k.SetBound_Symmetry(thread);
            w.SetBound_Symmetry(thread);
        }
        public void SetBound_Flow_ByTurbIntenAndVisRate(IField<double> TurbulenceIntensity, IField<double> TurbRadio, FvBound thread)
        {
            throw new NotImplementedException();
            //BoundkValue boundk = new BoundkValue();
            //boundk.SetParameter(TurbulenceIntensity, VelVector.VelocityMagnitude);
            //BoundwValue boundw = new BoundwValue();
            //boundw.SetParameter(TurbRadio, k_x, Density, DynamicViscosity);
            //SetBound_Flow(boundk, boundw, thread);
        }
        public void SetBound_Flow(IField<double> k_flowValue, IField<double> w_flowValue, FvBound thread)
        {
            k.SetBound_FlowValue(k_flowValue, thread);
            w.SetBound_FlowValue(w_flowValue, thread);
        }
        public void Cal_MomGamma(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //{
            //    (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            //    var den = Density.GetCellVars(Zone.Inner);
            //    for (; i1 < endp1; i1++)
            //    {
            //        //这个公式待确认，紧邻边界的网格也是这样算吗？！！！！！！！
            //        //边界处的速度梯度很大，计算出的粘度会很大
            //        MatrixR3C3 dudx = calDuDx(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1]);
            //        SymMatrixR3C3 Sij_ = calFavreAveragedStrainRate(dudx);//不可压，忽略后一项
            //        double wbo = Math.Max(w_x[i1], Clim * Math.Sqrt(two_div_betastar * Sij_.LengthSquared()));
            //        miu_T[i1] = den.GetCellVar(i1) * k_x[i1] / wbo;
            //    }
            //}
            //MultiThreads.Sync(threadInfo);
            //calGamma_byAverage(Zone.Inner, threadInfo);
            //foreach (var b in VelVector.NoSlipWalls)
            //{
            //    (int i1, int endp1) = b.InBound.FaceIdx.Split(threadInfo);
            //    double[] GammaCommon = TotalViscosity.CommonArray;
            //    for (; i1 < endp1; i1++)
            //    {
            //        GammaCommon[i1] = BoundValueCal.Cal_miuT(i1, ref Faces[i1]);
            //    }
            //}
            //foreach (var b in VelVector.Symmetries)
            //{
            //    calGamma_byBoundOwner(b.InBound, threadInfo);
            //}
            //foreach (var b in VelVector.VelocityInlets)
            //{
            //    calGamma_byBoundOwner(b.InBound, threadInfo);
            //}
            //foreach (var b in VelVector.PressureInOuts)
            //{
            //    calGamma_byBoundOwner(b.InBound, threadInfo);
            //}
        }

        public void Cal_kGamma(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //cal_kwVis(Zone, sigmastar, threadInfo);
            //MultiThreads.Sync(threadInfo);
            //calGamma_byAverage(Zone.Inner, threadInfo);
            //foreach (var b in VelVector.VelocityInlets)
            //{
            //    calGamma_byBoundNeighbor(b.InBound, threadInfo);
            //}
            //foreach (var b in VelVector.PressureInOuts)
            //{
            //    calGamma_byBoundNeighbor(b.InBound, threadInfo);
            //}
        }

        public void Add_k_Ab_byConstitutiveEquation(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //(int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);

            //var den = Density.GetCellVars(Zone.Inner);
            //bool[] isWallAdj = WallDistance.IsCellAdjWall;
            //double[] k_b = k.B;
            //double[] k_A_Diag = k.ADiag;
            //for (; i1 < endp1; i1++)
            //{
            //    ref FvCell cell = ref Cells[i1];
            //    double des = den.GetCellVar(i1);
            //    double w = w_x[i1];
            //    if (isWallAdj[i1]) //判断是否边界临近壁面
            //    {
            //        double Pk = double.MaxValue;
            //        for (int i2 = cell.InnerOFCount; i2 < cell.OFaces.Length; i2++)
            //        {//对于多个壁面的情况，取最小值？？
            //            int fid = cell.OFaces[i2];
            //            ref FvFace face = ref Faces[fid];
            //            if (isWallAdj[face.NSideCell])
            //            {
            //                double fPk = BoundValueCal.CalPk(fid, ref face);
            //                if (fPk < Pk) { Pk = fPk; }
            //            }
            //        }
            //        if (Pk == double.MaxValue) throw new Exception("湍流边界计算错误");
            //        k_b[i1] += Pk * cell.Volumn;
            //    }
            //    else
            //    {
            //        MatrixR3C3 dudx = calDuDx(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1]);
            //        SymMatrixR3C3 Sij_ = calFavreAveragedStrainRate(dudx);//不可压，忽略后一项
            //        SymMatrixR3C3 rhotau = 2 * miu_T[i1] * Sij_ - new DiagonalMatrixR3C3(LxTool.TwoThird * des * k_x[i1]);
            //        double Pk = MatrixR3C3.IJProduct(rhotau, dudx);
            //        k_b[i1] += Pk * cell.Volumn;
            //    }
            //    k_A_Diag[i1] += betastar * des * w * cell.Volumn;
            //}
        }
        protected MatrixR3C3 calDuDx(ref Vector3d gradU, ref Vector3d gradV, ref Vector3d gradW)
        {
            return new MatrixR3C3(
                gradU.X, gradU.Y, gradU.Z,
                gradV.X, gradV.Y, gradV.Z,
                gradW.X, gradW.Y, gradW.Z);
        }
        protected SymMatrixR3C3 calFavreAveragedStrainRate(in MatrixR3C3 dudx)//对于WilCox论文Sij
        {
            return new SymMatrixR3C3(dudx.a11,
                0.5 * (dudx.a21 + dudx.a12), dudx.a22,
                0.5 * (dudx.a31 + dudx.a13), 0.5 * (dudx.a32 + dudx.a23), dudx.a33);
        }
        protected SymMatrixR3C3 calZeroTraceFavreAveragedStrainRate(in MatrixR3C3 dudx)//对于WilCox论文Sij带横杠
        {
            double temp = (dudx.a11 + dudx.a22 + dudx.a33) / 3;
            return new SymMatrixR3C3(dudx.a11 - temp,
                0.5 * (dudx.a21 + dudx.a12), dudx.a22 - temp,
                0.5 * (dudx.a31 + dudx.a13), 0.5 * (dudx.a32 + dudx.a23), dudx.a33 - temp);
        }
        protected SymMatrixR3C3 calGalileanFavreAveragedStrainRate(in MatrixR3C3 dudx)//对于WilCox论文Sij带尖
        {
            double temp = (dudx.a11 + dudx.a22 + dudx.a33) / 2;
            return new SymMatrixR3C3(dudx.a11 - temp,
                0.5 * (dudx.a21 + dudx.a12), dudx.a22 - temp,
                0.5 * (dudx.a31 + dudx.a13), 0.5 * (dudx.a32 + dudx.a23), dudx.a33 - temp);
        }
        public void Cal_wGamma(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //cal_kwVis(Zone, sigma, threadInfo);
            //MultiThreads.Sync(threadInfo);
            //calGamma_byAverage(Zone.Inner, threadInfo);
            //foreach (var b in VelVector.VelocityInlets)
            //{
            //    calGamma_byBoundNeighbor(b.InBound, threadInfo);
            //}
            //foreach (var b in VelVector.PressureInOuts)
            //{
            //    calGamma_byBoundNeighbor(b.InBound, threadInfo);
            //}
        }
        public void Set_w_Ab_byConstitutiveEquation_AfterSetAll(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //(int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);
            //var den = Density.GetCellVars(Zone.Inner);
            //bool[] isWallAdj = WallDistance.IsCellAdjWall;
            //double[] w_b = w.B;
            //double[] w_A_Diag = w.ADiag;
            //Vector3d[] w_grad = w.Grad;
            //double[] w_A_LU_O = w.ALUO;
            //double[] w_A_LU_N = w.ALUN;
            //Vector3d[] k_grad = k.Grad;
            //for (; i1 < endp1; i1++)
            //{
            //    ref FvCell cell = ref Cells[i1];
            //    double des = den.GetCellVar(i1);
            //    double k = k_x[i1];
            //    double w = w_x[i1];
            //    if (isWallAdj[i1])
            //    {
            //        double omegaC = 0;
            //        for (int i2 = cell.InnerOFCount; i2 < cell.OFaces.Length; i2++)
            //        {//对于多个壁面的情况，取最大值？？
            //            int fid = cell.OFaces[i2];
            //            ref FvFace face = ref Faces[fid];
            //            if (isWallAdj[face.NSideCell])
            //            {
            //                double fw = BoundValueCal.Cal_omegaC(fid, ref face);
            //                if (fw > omegaC) { omegaC = fw; }
            //            }
            //        }
            //        if (omegaC == 0) throw new Exception("湍流边界计算错误");
            //        w_b[i1] = omegaC * w_A_Diag[i1];
            //        w_x[i1] = omegaC;
            //        foreach (int f in cell.OFaces)
            //        {
            //            w_A_LU_O[f] = 0;
            //        }
            //        foreach (int f in cell.NFaces)
            //        {
            //            w_A_LU_N[f] = 0;
            //        }
            //    }
            //    else
            //    {
            //        MatrixR3C3 dudx = calDuDx(ref U_grad[i1], ref V_grad[i1], ref W_grad[i1]);
            //        SymMatrixR3C3 Sij_ = calFavreAveragedStrainRate(dudx);//不可压，忽略后一项
            //        SymMatrixR3C3 rhotau = 2 * miu_T[i1] * Sij_ - new DiagonalMatrixR3C3(LxTool.TwoThird * des * k_x[i1]);
            //        double Pk = MatrixR3C3.IJProduct(rhotau, dudx);//到这里与k方程中的源项相同
            //        double dkdw = Vector3d.DotProduct(k_grad[i1], w_grad[i1]);
            //        double sigmaddkdw = dkdw > 0 ? sigmado * des / w : 0;
            //        w_b[i1] += (alpha * w / k * Pk + sigmaddkdw) * cell.Volumn;
            //        double OOS = calOmegaOmegaSijjian(dudx, Sij_);
            //        double Chiw = Math.Abs(OOS / Math.Pow(betastar * w, 3));
            //        double beta = betao * ((1 + 85 * Chiw) / (1 + 100 * Chiw));
            //        w_A_Diag[i1] += beta * des * w * cell.Volumn;
            //    }
            //}
        }

        protected void cal_kwVis(IIndexHost indexHost, double sigmaStarForkAndsigmaForw, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = indexHost.CellIdx.Split(threadInfo);
            var den = Density.GetCellVarArr(indexHost);
            for (; i1 < endp1; i1++)
            {
                miu_T[i1] = sigmaStarForkAndsigmaForw * den.GetCellVar(i1) * k_x[i1] / w_x[i1];
            }
        }

        protected void calGamma_byAverage(IIndexHost indexHost, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] GammaCommon = TotalViscosity.CommonArray;
            //(int i1, int endp1) = indexHost.GetFaceIdx().Split(threadInfo);
            //var dyn = DynamicViscosity.GetFaceVars(indexHost);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dyn.GetFaceVar(i1) + Faces[i1].CalAverage(miu_T);
            //}
        }
        protected void calGamma_byBoundOwner(FvBound bound, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] GammaCommon = TotalViscosity.CommonArray;
            //(int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            //var dyn = DynamicViscosity.GetFaceVars(bound);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dyn.GetFaceVar(i1) + miu_T[Faces[i1].OSideCell];
            //}
        }
        protected void calGamma_byBoundNeighbor(FvBound bound, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] GammaCommon = TotalViscosity.CommonArray;
            //(int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            //var dyn = DynamicViscosity.GetFaceVars(bound);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = dyn.GetFaceVar(i1) + miu_T[Faces[i1].NSideCell];
            //}
        }
        protected double calOmegaOmegaSijjian(in MatrixR3C3 dudx, in SymMatrixR3C3 Sijjian)
        {
            double a21 = 0.5 * (dudx.a21 - dudx.a12);
            double a31 = 0.5 * (dudx.a31 - dudx.a13);
            double a32 = 0.5 * (dudx.a32 - dudx.a23);
            return -a21 * a21 * (Sijjian.a11 + Sijjian.a22)
                - a31 * a31 * (Sijjian.a11 + Sijjian.a33)
                - a32 * a32 * (Sijjian.a33 + Sijjian.a22)
                + 2 * (a21 * a32 * Sijjian.a31 - a31 * a32 * Sijjian.a21 - a21 * a31 * Sijjian.a32);
        }
        public abstract class CalBoundBase
        {
            protected IField<double> density;
            protected IField<double> dynVis;
            protected SIMPLEVelocity VelVector;
            protected double[] WallDistance;
            protected double[] k_x;
            public IField<double> ThermalConductivity;
            public IField<double> SpecificHeat;
            public void SetTurb(TurbulencekOmega inTurb)
            {
                this.VelVector = inTurb.VelVector;
                this.density = inTurb.Density;
                this.WallDistance = inTurb.WallDistance.InnDis;
                this.dynVis = inTurb.DynamicViscosity;
                this.k_x = inTurb.k.InnX;
            }
            public void SetThermal(IField<double> ThermalConductivity, IField<double> SpecificHeat)
            {
                this.ThermalConductivity = ThermalConductivity;
                this.SpecificHeat = SpecificHeat;
            }
            protected double cal_u(ref FvInnFace face)
            {
                throw new NotImplementedException();
                //Vector3d vel = VelVector.GetVelocity(face.OSideCell) - VelVector.GetVelocity(face.NSideCell);
                //Vector3d vel_parallel = Vector3d.ProjectVectorToPlane(vel, face.Normal);
                //return vel_parallel.Length();
            }
            public virtual double Cal_miuT(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
            }
            public virtual double CalPk(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
            }
            public virtual double Cal_omegaC(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
            }
            public virtual double Cal_ThermalConductivity(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
            }
        }
        public class LowYPlusBound : CalBoundBase
        {
            //public override double Cal_miuT(int fid, ref FvFace face)
            //{
            //    return dynVis.GetFaceValue(fid);
            //}
            //public override double CalPk(int fid, ref FvFace face)
            //{
            //    double dudy = cal_u(ref face) / WallDistance[face.OSideCell];
            //    return dynVis.GetFaceValue(fid) * dudy * dudy;
            //}
            //public override double Cal_omegaC(int fid, ref FvFace face)
            //{
            //    double y = WallDistance[face.OSideCell];
            //    return 6 * dynVis.GetFaceValue(fid) / (Cbeta1 * density.GetCellValue(face.OSideCell) * y * y);
            //}
            //public override double Cal_ThermalConductivity(int fid, ref FvFace face)
            //{
            //    return ThermalConductivity.GetFaceValue(fid);
            //}
        }
        public class HighYPlusBound : CalBoundBase
        {
            protected double cal_ustar(out double kSqrted, double k)
            {
                kSqrted = Math.Sqrt(k);
                return Cmu14 * kSqrted;
            }
            public override double Cal_miuT(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
                //double k = k_x[face.OSideCell];
                //double kSqrted = Math.Sqrt(k);
                //double uStar = Cmu14 * kSqrted;
                //double y = WallDistance[face.OSideCell];
                //double rho = density.GetCellValue(face.OSideCell);
                //double yPlus = y * rho * uStar / dynVis.GetFaceValue(fid);
                //yPlus = Math.Max(yPlus, 11);//yplus太小时计算ln函数会出错
                //return rho * uStar * y / (Karman_Reciprocal * Math.Log(yPlus) + B);
            }
            public override double CalPk(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
                //double rho = density.GetCellValue(face.OSideCell);
                //double k = k_x[face.OSideCell];
                //double kSqrted = Math.Sqrt(k);
                //double uStar = Cmu14 * kSqrted;
                //double y = WallDistance[face.OSideCell];
                //return rho * Math.Pow(uStar, 3) / (Karman * y);
            }
            public override double Cal_omegaC(int fid, ref FvInnFace face)
            {
                double k = k_x[face.OSideCell];
                double kSqrted = Math.Sqrt(k);
                double y = WallDistance[face.OSideCell];
                return kSqrted / (Karman * Cmu14 * y);
            }
            public override double Cal_ThermalConductivity(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
                //double miu = dynVis.GetFaceValue(fid);
                //double k = k_x[face.OSideCell];
                //double kSqrted = Math.Sqrt(k);
                //double uStar = Cmu14 * kSqrted;
                //double y = WallDistance[face.OSideCell];
                //double rho = density.GetCellValue(face.OSideCell);
                //double yPlus = y * rho * uStar / miu;
                //yPlus = Math.Max(yPlus, 11);//yplus太小时计算ln函数会出错
                //double specificHeat = SpecificHeat.GetCellValue(face.OSideCell);
                //double Pr = specificHeat * miu / ThermalConductivity.GetCellValue(face.OSideCell);
                //double Pr_div_Prt = Pr * one_div_PrT;
                //double P = 9.24 * (Math.Pow(Pr_div_Prt, 0.75) - 1) * (1 + 0.28 * Math.Exp(-0.007 * Pr_div_Prt));
                //double TPlus = PrT * (Karman_Reciprocal * Math.Log(yPlus) + P);
                //return rho * specificHeat * uStar * y / TPlus;
            }
        }
        public class EnhanceYPlusBound : CalBoundBase
        {
            protected void calyrhomu(out double y, out double rho, out double mu, out double yrho_miu, int faceIndex, ref FvInnFace face)
            {
                throw new NotImplementedException();
                //rho = density.GetCellValue(face.OSideCell);
                //mu = dynVis.GetFaceValue(faceIndex);
                //y = WallDistance[face.OSideCell];//第一层cell
                //yrho_miu = y * rho / mu;
            }
            protected void calRedgamma(out double gamma, out double kSqrted, ref FvInnFace face, double yrho_mu)
            {
                kSqrted = Math.Sqrt(k_x[face.OSideCell]);
                double Red = kSqrted * yrho_mu;//这个地方用的Face上的Gamma，为了统一，Gamma都是用的Face
                gamma = Math.Exp(-Red / 11.0);
            }
            protected void calustar(out double uStar, out double u, out double yPlus, ref FvInnFace face, double gamma, double yrho_mu, double kSqrted)
            {
                u = cal_u(ref face);
                uStar = gamma * Math.Sqrt(u / yrho_mu) + (1 - gamma) * Cmu14 * kSqrted;//感觉应该只包括后一项
                yPlus = yrho_mu * uStar;
            }
            protected double calBound_mu_w(int faceIndex, ref FvInnFace face)
            {
                calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
                calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
                calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
                //double uPlus = gamma * yPlus + gamma_1 * (Karman_Reciprocal * Math.Log(yPlus) + B);//必须采用更详细的公式，这样弄第二项可能为负//也用了混合函数，StarCCM+帮助中，由更详细的公式，后续可参考改进
                double uPlus = yPlus < 11 ? yPlus : Karman_Reciprocal * Math.Log(yPlus) + B;
                return gamma * mu + (1 - gamma) * (rho * uStar * y) / uPlus;//壁面处粘性系数mu
            }
            protected double calBound_Pk(int faceIndex, ref FvInnFace face)
            {
                calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
                calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
                calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
                return gamma * mu * Math.Pow(u / y, 2) + (1 - gamma) * rho * Math.Pow(uStar, 3) / (Karman * y);
            }
            protected double calBound_omegaC(int faceIndex, ref FvInnFace face)
            {
                calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, faceIndex, ref face);
                calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
                return gamma * 60 * mu / (Cbeta1 * rho * Math.Pow(y, 2)) + (1 - gamma) * (kSqrted / (Karman * Cmu14 * y));
            }
            public override double Cal_ThermalConductivity(int fid, ref FvInnFace face)
            {
                throw new NotImplementedException();
                //calyrhomu(out double y, out double rho, out double mu, out double yrho_mu, fid, ref face);
                //calRedgamma(out double gamma, out double kSqrted, ref face, yrho_mu);
                //calustar(out double uStar, out double u, out double yPlus, ref face, gamma, yrho_mu, kSqrted);
                //double specificHeat = SpecificHeat.GetCellValue(face.OSideCell);
                //double Pr = specificHeat * mu / ThermalConductivity.GetCellValue(face.OSideCell);
                //double Gamma = 0.01 * Math.Pow(Pr * yPlus, 4) / (1 + 5 * Math.Pow(Pr, 3) * yPlus);
                //double exp_Gamma = Math.Exp(-Gamma);
                //double Pr_div_Prt = Pr * one_div_PrT;
                //double P = 9.24 * (Math.Pow(Pr_div_Prt, 0.75) - 1) * (1 + 0.28 * Math.Exp(-0.007 * Pr_div_Prt));
                //double TPlus = exp_Gamma * Pr * yPlus + (1 - exp_Gamma) * PrT * (Karman_Reciprocal * Math.Log(yPlus) + P);
                //return rho * specificHeat * uStar * y / TPlus;
            }
        }

        public class BoundkValue : FieldBase<double>
        {
            public IField<double> TurbulenceIntensity;// TurbIntensitySquared_div_2;
            protected IField<double> Velocity;
            public void SetParameter(IField<double> turbIntensity, IField<double> VelVector)
            {
                // this.TurbIntensitySquared_div_2 = 0.5 * Math.Pow(turbIntensity, 2);
                this.TurbulenceIntensity = turbIntensity;
                this.Velocity = VelVector;
            }
            //public double GetCellValue(int cellIndex)
            //{
            //    double I = TurbulenceIntensity.GetCellValue(cellIndex);
            //    double v = Velocity.GetCellValue(cellIndex);
            //    return 0.5 * (I * I) * (v * v);
            //}
            public override ICellVarArr<double> GetCellVarArr(IIndexHost indexHost)
            {
                return new Boundk()
                {
                    TurbulenceIntensity = TurbulenceIntensity.GetCellVarArr(indexHost),
                    Velocity = Velocity.GetCellVarArr(indexHost)
                };
            }
            public struct Boundk : ICellVarArr<double>
            {
                public ICellVarArr<double> TurbulenceIntensity;// TurbIntensitySquared_div_2;
                public ICellVarArr<double> Velocity;
                public double GetCellVar(int cellID)
                {
                    double I = TurbulenceIntensity.GetCellVar(cellID);
                    double v = Velocity.GetCellVar(cellID);
                    return 0.5 * (I * I) * (v * v);
                }
            }
        }
        public class BoundwValue : FieldBase<double>
        {
            IField<double> VisRatio;
            double[] k_x;
            IField<double> density;
            IField<double> dynVis;
            public void SetParameter(IField<double> VisRatio, double[] k_x, IField<double> density, IField<double> dynVis)
            {
                this.VisRatio = VisRatio;
                this.k_x = k_x;
                this.density = density;
                this.dynVis = dynVis;
            }
            public override ICellVarArr<double> GetCellVarArr(IIndexHost indexHost)
            {
                return new Boundw()
                {
                    VisRatio = VisRatio.GetCellVarArr(indexHost),
                    k_x = k_x,
                    density = density.GetCellVarArr(indexHost),
                    dynVis = dynVis.GetCellVarArr(indexHost),
                };
            }
            public struct Boundw : ICellVarArr<double>
            {
                public ICellVarArr<double> VisRatio;
                public double[] k_x;
                public ICellVarArr<double> density;
                public ICellVarArr<double> dynVis;
                public double GetCellVar(int cellIndex)
                {
                    return density.GetCellVar(cellIndex) * k_x[cellIndex] / dynVis.GetCellVar(cellIndex) * VisRatio.GetCellVar(cellIndex);
                }
            }
        }

        public IField<double> ThermalConductivity;
        public IField<double> SpecificHeat;
        public void SetZone_part3_TemperatureGamma(bool isEnableTemperature, IField<double> specificHeat, IField<double> thermalConductivity)
        {
            if (!isEnableTemperature)
            {
                this.SpecificHeat = null;
                this.ThermalConductivity = null;
                BoundValueCal.SetThermal(null, null);
                return;
            }
            this.SpecificHeat = specificHeat;
            this.ThermalConductivity = thermalConductivity;
            BoundValueCal.SetThermal(thermalConductivity, specificHeat);
        }
        public void Cal_TemperatureGamma(ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] GammaCommon = TotalViscosity.CommonArray;
            //{
            //    InnerIndex.FaceIDs.GetSeg(  threadInfo);
            //    for (; i1 < endp1; i1++)
            //    {
            //        GammaCommon[i1] = ThermalConductivity.GetFaceValue(i1)
            //            + SpecificHeat.GetFaceValue(i1) * one_div_PrT * Faces[i1].CalAverage(miu_T);
            //    }
            //}
            //foreach (var b in VelVector.NoSlipWalls)
            //{
            //    b.BoundIndex.FaceIDs.GetSeg(  threadInfo);
            //    for (; i1 < endp1; i1++)
            //    {
            //        GammaCommon[i1] = BoundValueCal.Cal_ThermalConductivity(i1, ref Faces[i1]);
            //    }
            //}
            //foreach (var b in VelVector.Symmetries)
            //{
            //    b.BoundIndex.FaceIDs.GetSeg(  threadInfo);
            //    Array.Clear(GammaCommon, i1, endp1 - i1);
            //}
            //foreach (var b in VelVector.VelocityInlets)
            //{
            //    calBound_ThermalConductivity(b.BoundIndex, threadInfo);
            //}
            //foreach (var b in VelVector.PressureInOuts)
            //{
            //    calBound_ThermalConductivity(b.BoundIndex, threadInfo);
            //}
        }
        protected void calBound_ThermalConductivity(FvBound bound, ThreadInfo threadInfo)
        {
            throw new NotImplementedException();
            //double[] GammaCommon = TotalViscosity.CommonArray;
            //indexNumber.FaceIDs.GetSeg(  threadInfo);
            //for (; i1 < endp1; i1++)
            //{
            //    GammaCommon[i1] = ThermalConductivity.GetFaceValue(i1)
            //        + SpecificHeat.GetFaceValue(i1) * one_div_PrT * miu_T[Faces[i1].OSideCell];
            //}
        }

    }
}
